Event History Analysis II — Cox proportional hazard regression
Dagen program
Øvelsen til forelæsning 3.
Data og data management
Cox-regression
Analyse i praksis
Øvelse
Forventet udbytte
Viden
Cox-regressionen og dens forudsætninger.
Cox-regressionens fordele (og ulemper).
Færdigheder
Kodning af tidsvarierende variable til tidskonstante variable.
Grafisk præsentation af resultater fra en cox-regression.
Kompetancer
Analyse af Cox-regression i R.
Subject-Period file og single-episode file
Hvor vi sidste gang arbejde med en single-episode file, arbejder vi i konteksten af regressionsanalyse med mikrodata, ofte med en Subject-Period file — dvs. en fil, hvor hvert individ fremgår flere gange over tid.
Konstruktion af event file
Den centrale kode-mæssige udfordring er at konstruere event variablen.
Funktionerne lag() og lead() er centrale, når vi skal konstruere vores event variabel. Disse funktioner bruges til at hente data fra henholdsvis den foregående og den efterfølgende (observations)række i vores dataframe. Altså, lag() og lead() bruges til at “kalde” oplysninger fra henholdsvis før og efter en given observationsrække.
Spg.: Hvordan løser det vores “problem” med at skabe vores event indikator?
02:00
Indikator forstartstidspunktetper person:
library(tidyverse)df.rå <-read_excel("eksempel_data.xlsx")df.lagged <- df.rå %>%mutate(person_start =case_when(lag(id) == id ~0, TRUE~1))
Konstruér starttidspunkt i forløbet
Fiktive data
id
år
børn
alder
region18
tid
event
person_start
1
1
2000
0
18
Nord
1
0
1
1
2001
0
19
Nord
2
0
0
1
2002
1
20
Syd
3
1
0
1
2003
1
21
Syd
4
0
0
1
2004
1
22
Syd
5
0
0
2
2
2003
0
20
Syd
3
0
1
2
2004
0
21
Øst
4
0
0
2
2005
0
22
Øst
5
0
0
2
2006
0
23
Nord
6
0
0
2
2007
0
24
Nord
7
0
0
3
3
1999
0
18
Nord
1
0
1
3
2000
0
19
Nord
2
0
0
3
2001
0
20
Vest
3
0
0
3
2002
0
21
Vest
4
0
0
3
2003
1
22
Nord
5
1
0
4
4
1993
0
23
Øst
6
0
1
4
1994
0
24
Øst
7
0
0
4
1995
0
25
Øst
8
0
0
4
1996
1
26
Øst
9
1
0
4
1997
1
27
Øst
10
0
0
Indikator forsluttidspunktetper person:
library(tidyverse)df.rå <-read_excel("eksempel_data.xlsx")df.lagged <- df.rå %>%mutate(person_start =case_when(lag(id) == id ~0, TRUE~1)) %>%mutate(person_slut =case_when(lead(id) == id ~0, TRUE~1))
Konstruér starttidspunkt i forløbet
Fiktive data
id
år
børn
alder
region18
tid
event
person_start
person_slut
1
1
2000
0
18
Nord
1
0
1
0
1
2001
0
19
Nord
2
0
0
0
1
2002
1
20
Syd
3
1
0
0
1
2003
1
21
Syd
4
0
0
0
1
2004
1
22
Syd
5
0
0
1
2
2
2003
0
20
Syd
3
0
1
0
2
2004
0
21
Øst
4
0
0
0
2
2005
0
22
Øst
5
0
0
0
2
2006
0
23
Nord
6
0
0
0
2
2007
0
24
Nord
7
0
0
1
3
3
1999
0
18
Nord
1
0
1
0
3
2000
0
19
Nord
2
0
0
0
3
2001
0
20
Vest
3
0
0
0
3
2002
0
21
Vest
4
0
0
0
3
2003
1
22
Nord
5
1
0
1
4
4
1993
0
23
Øst
6
0
1
0
4
1994
0
24
Øst
7
0
0
0
4
1995
0
25
Øst
8
0
0
0
4
1996
1
26
Øst
9
1
0
0
4
1997
1
27
Øst
10
0
0
1
Hvis vi nu vil definere begivenheden “fødsel af første barn” er udfordringen at “lokalisere” der hvor kvinden går fra 0 børn til 1 barn.
Jeg bruger her både case_when() og if_else() for illustrationen til begge funktinoer.
library(tidyverse)df.rå <-read_excel("eksempel_data.xlsx")df.lagged <- df.rå %>%mutate(person_start =case_when(lag(id) == id ~0, TRUE~1)) %>%mutate(person_slut =case_when(lead(id) == id ~0, TRUE~1)) %>%group_by() %>%mutate(begivenhed =if_else(lag(børn) ==0& børn ==1& person_start ==0, 1, 0)) %>%ungroup()
Konstruér indikator for begivenhed
Fiktive data
id
år
børn
alder
region18
tid
event
person_start
person_slut
begivenhed
1
1
2000
0
18
Nord
1
0
1
0
0
1
2001
0
19
Nord
2
0
0
0
0
1
2002
1
20
Syd
3
1
0
0
1
1
2003
1
21
Syd
4
0
0
0
0
1
2004
1
22
Syd
5
0
0
1
0
2
2
2003
0
20
Syd
3
0
1
0
0
2
2004
0
21
Øst
4
0
0
0
0
2
2005
0
22
Øst
5
0
0
0
0
2
2006
0
23
Nord
6
0
0
0
0
2
2007
0
24
Nord
7
0
0
1
0
3
3
1999
0
18
Nord
1
0
1
0
0
3
2000
0
19
Nord
2
0
0
0
0
3
2001
0
20
Vest
3
0
0
0
0
3
2002
0
21
Vest
4
0
0
0
0
3
2003
1
22
Nord
5
1
0
1
1
4
4
1993
0
23
Øst
6
0
1
0
0
4
1994
0
24
Øst
7
0
0
0
0
4
1995
0
25
Øst
8
0
0
0
0
4
1996
1
26
Øst
9
1
0
0
1
4
1997
1
27
Øst
10
0
0
1
0
library(tidyverse)df.rå <-read_excel("eksempel_data.xlsx")df.lagged <- df.rå %>%mutate(person_start =case_when(lag(id) == id ~0, TRUE~1),person_slut =case_when(lead(id) == id ~0, TRUE~1)) %>%group_by(id) %>%mutate(begivenhed =if_else(lag(børn) ==0& børn ==1& person_start ==0, 1, 0),periode_slut =case_when(begivenhed ==1| børn ==0& person_slut ==1~1, TRUE~0)) %>%slice(1:which.max(periode_slut ==1)) %>%ungroup()
Konstruér subject-period file
Fiktive data
id
år
børn
alder
region18
tid
event
person_start
person_slut
begivenhed
periode_slut
1
1
2000
0
18
Nord
1
0
1
0
0
0
1
2001
0
19
Nord
2
0
0
0
0
0
1
2002
1
20
Syd
3
1
0
0
1
1
2
2
2003
0
20
Syd
3
0
1
0
0
0
2
2004
0
21
Øst
4
0
0
0
0
0
2
2005
0
22
Øst
5
0
0
0
0
0
2
2006
0
23
Nord
6
0
0
0
0
0
2
2007
0
24
Nord
7
0
0
1
0
1
3
3
1999
0
18
Nord
1
0
1
0
0
0
3
2000
0
19
Nord
2
0
0
0
0
0
3
2001
0
20
Vest
3
0
0
0
0
0
3
2002
0
21
Vest
4
0
0
0
0
0
3
2003
1
22
Nord
5
1
0
1
1
1
4
4
1993
0
23
Øst
6
0
1
0
0
0
4
1994
0
24
Øst
7
0
0
0
0
0
4
1995
0
25
Øst
8
0
0
0
0
0
4
1996
1
26
Øst
9
1
0
0
1
1
Herfra er en Single-episode file ligefrem at konstruere, ved:
library(tidyverse)df.rå <-read_excel("eksempel_data.xlsx")df.lagged <- df.rå %>%mutate(person_start =case_when(lag(id) == id ~0, TRUE~1)) %>%mutate(person_slut =case_when(lead(id) == id ~0, TRUE~1)) %>%group_by(id) %>%mutate(begivenhed =if_else(lag(børn) ==0& børn ==1& person_start ==0, 1, 0)) %>%mutate(periode_slut =case_when(begivenhed ==1| børn ==0& person_slut ==1~1, TRUE~0)) %>%ungroup() %>%filter(periode_slut ==1)
Konstruér single-episode file
Fiktive data
grp
id
år
børn
alder
region18
tid
event
person_start
person_slut
begivenhed
periode_slut
1
1
2002
1
20
Syd
3
1
0
0
1
1
2
2
2007
0
24
Nord
7
0
0
1
0
1
3
3
2003
1
22
Nord
5
1
0
1
1
1
4
4
1996
1
26
Øst
9
1
0
0
1
1
Konstruktion af event file
Hvad så med en Single-episode file i stil med First_child.rda. Hvad er fordele og begrænsninger ved denne?
Cox PH-modellen er semi-parametrisk, hvilket betyder, at vi ikke behøver at specificere den underliggende sandsynlighedstæthedsfunktion (pdf) for hændelsestiden. (Genbesøg noter og slides til lektion 2, for introduktion til relevante pdf.) Dette gør modellen fleksibel, da vi kun estimerer effekten af kovariaterne på hazardraten, mens baseline hazard forbliver ukendt.
Modellen er robust og kan anvendes til data. Uanset om den underliggende fordeling for hændelsestiderne er Weibull, eksponentiel, lognormal eller en anden form, vil Cox-modellen give gode approksimationer, så længe den proportionale hazard-antagelse er opfyldt.
En væsentlig forudsætning for modellen er, at hazardraten for de forskellige grupper er proportionale over tid. Denne antagelse skal testes, da overtrædelser kan påvirke modellens anvendelighed.
En af de primøre fordele ved Cox-modellen er, at den nemt kan udvides til at inkludere tidsvarierende kovariater.
Cox proportional hazard model med fixed covariater.
For den \(i\)’te person til tiden \(t\) specificeres hazardfunktionen som:
\(h_{0}(t)\) (baseline hazard) angiver den underliggende risiko ved tid tt for en referenceperson, dvs. for en person med \(x_{1,i} = x_{2,i} = \dots = x_{k,i} = 0\). \(h_{0}(t)\) specificeres ikke a priori, hvilket betyder, at vi ikke antager nogen konkret form for dens fordeling.
\(\exp(\beta_{1} x_{1} + \dots + \beta_{k} x_{k})\) er den multiplicative effekt af de \(k\) covariater på hazardraten.
Bemærk at der ikke er noget egentlig intercept, da den faste effekt, som et intercept normalt ville repræsentere, er integreret i \(h_{0}(t)\). Ved estimering (via partial likelihood) udtrækkes koefficienterne \(\beta\), mens \(h_{0}(t)\) forbliver uspecificeret.
Cox proportional hazard model med fixed covariater.
For den \(i\)’te person til tiden \(t\) specificeres hazardfunktionen som:
Cox proportional hazard model med fixed covariater. Sammenligning af grupper.
Med afsæt i data fra sidste øvelse, First_child.rda, antager vi, at vi har en gruppeindikator \(x\), som defineres som en dummyvariabel:
\[
x =
\left\{\begin{matrix}
0, \: \text{hvis personen tilhører referencegruppen} \:\:\:\:\:\:\:\:
\\
1, \: \text{hvis personen ikke tilhører referencegruppen}
\end{matrix}\right.
\]
Konceptuelt svarer HR til Odds Ratio (OR) i logistisk regression, idet den udtrykker den relative ændring i hazardraten mellem grupper:
Hvis HR er \(=1\) har personen uden for referencegruppen samme hazard-rate som personen i referencegruppen.
Hvis HR er \(<1\) har personen uden for referencegruppen en lavere hazard-rate end personen i referencegruppen.
Hvis HR er \(>1\) har ersonen uden for referencegruppen en højere hazard-rate end personen i referencegruppen.
Den proportionale hazard model
Fordi \(h_{0}(t)\) bliver ophævet i tæller og nævner, gælder det at hazardrater mellem to personer ved samme tid er uafhængig af selve tiden. Det skyldes, at baseline hazardfunktionen \(h_{0}(t)\) ophæves, når vi danner forholdet (hazard-ratioen) mellem to individer. Dermed er hazard-ratioen udelukkende bestemt af de eksponentierede kovariat-effekter og antages at være konstant over tid.
Den proportionale hazard model
Fordi hazard-ratioen (HR) mellem to personer er konstant over tid, gælder det at:
\[
\hat{HR} =
\frac{h_{i}(t)}{h_{j}(t)} = \theta
\] hvor \(\theta\) er en konstant, der udelukkende bestemmes af de eksponentierede kovariateffekter.
Den proportionale hazard model
Beskrevet andeledes: alle starter med den samme “basisrisiko” for at blive ramt af en begivenhed. Når vi sammenligner to personer, ser vi kun på, hvordan deres individuelle egenskaber gør, at deres risikoer afviger fra denne “basisrisiko” – og det forhold (hazard-ratioen) er konstant over tid.
Fx: Gruppe 1 har en højere risiko; og denne risiko er\(n\)højere end gruppe 2 over hele forløbet.
03:00
adjusted overlevelseskurver
Den centrale interesse er fortsat på at estimere overlevelseskurven, \(S(t)\). Som I ved fra de forrige lektioner er der en relation mellem \(S(t)\) og \(h(t)\) (genlæs noter eller slide til lektion 2+3).
Når vi inkluderer individuelle karakteristika (covariater) i analysen, udtrykker vi hazardfunktionen for en person med covariatvektoren \(\mathbf{X}\) som: \[
h(t, \mathbf{X}) = h_{j}(t) \times \text{exp}(\mathbf{X}\beta)
\]
hvor
\(h_{j}(t)\) er baseline hazard: den risiko, der gælder for referencegruppen
\(\text{exp}(\mathbf{X}\beta)\) beskriver, hvordan de enkelte covariater ændrer denne risiko.
adjusted overlevelseskurver
Dette giver os en “kontrolleret” overlevelsesfunktion, som justeres for covariaterne:
\(\hat{S}(t)\) er den “kontrollerede” \(S(t)\), der bestemmer overlevelsesforløbet for en person med specifikke karakteristika, hvor effekten af disse karakteristika er “kontrolleret” for.
library(tidyverse)library(survival)library(ggsurvfit)df <-read_csv(here::here("flytte.csv"))km.m1 <-survfit(Surv(dur, status) ~1, type ="kaplan-meier", data = df)km.m1 %>%ggsurvfit(color ="red") +add_confidence_interval()
Cox i praksis (sammenlignet med KM)
library(tidyverse)library(survival)library(ggsurvfit)df <-read_csv(here::here("flytte.csv"))km.m1 <-survfit(Surv(dur, status) ~ udd, type ="kaplan-meier", data = df)km.m1 %>%ggsurvfit() +add_confidence_interval()
Cox i praksis (sammenlignet med KM)
library(tidyverse)library(survival)library(ggsurvfit)df <-read_csv(here::here("flytte.csv"))km.m1 <-survfit(Surv(dur, status) ~ udd + aargang, type ="kaplan-meier", data = df)km.m1 %>%ggsurvfit() +add_confidence_interval()
Cox i praksis
Hvis vi vil vide hvordan overlevelse varierer i henhold til en specifik covariat, kan vi stratificere analysen og stadig “kontrollere” for baggrundsvariable. \(S(t)\) estimeres for hvert strata, som her er uddannelsesretning, kontrolleret for årgang.
Hvor bekendte er i med at “joine” data? Skal vi køre en introduktion inden øvelserne?
Dagens øvelse
Vi kigger nærmere på alder for fødsel af første barn (inklusiv alderen på faren ved fødsel).
load(".../forlob_06032023.rda")
Datasættet minder om øvelsesdata fra sidste gang.
Det adskiller sig ved at det
er person-period data (indeholder alle obs. frem til hændelse eller datastop)
hændelsesvariable skal dannes ud fra C_ANTBOERNF
regionsvariablen er nu tidsvaerierende
Datamanagement
Lav region (tidsvarierende variabel) om til en tidsinvariant variabel, der indikerer regionsbopælen i personens 18. leveår (forløbets start). Den skal ligne region18 fra sidste lektion. hint:*_join()
Konstruer ud fra C_ANTBOERNF en event-variabel: Tag højde for hvordan C_ANTBOERNF er defineret af DST (læs dokumentation på DST Times… Google)
Diskuter hvordan defininationen kan være et problem og hvordan vi kan løse problemt… hvis det er et problem?
Dan en variabel, der er en indentikator for den sidste observation i forløbet. Den kan enten være en hændelse eller en censurering.
Reducér data til kun at indeholde den sidste observation iforløbet.
Analyse
Lav en KM
Forskellige tests - analysér og diskutér
Lav et KM plot
Cox-regression (bivariat)
Bivariat med køn - analysér og diskutér hvad resultatet betyder
Cox-regression (multivariat)
Indsæt køn, region (tidsinvariant) og uddannelse
Fortolk resultaterne
Plot kønnenes overlevelse med udgangspunkt i cox-regressionen